Sandpile model on an optimized scale-free network on Euclidean space 
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Deterministic sandpile models are studied on a cost optimized Barabasi-Albert (BA) scale-free 
network whose nodes are the sites of a square lattice. For the optimized BA network, the sandpile 
model has the same critical behaviour as the BTW sandpile, whereas for the un-optimized BA 
network the critical behaviour is mean-field like. 
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Various models of Statistical Mechanics which are usu- 
ally studied on regular lattices are being studied in recent 
years on graphs or networks of very complex structures. 
For example, the Ising model has been studied both on 
the Small- World Networks (SWN) Q, and on Scale- 
Free Networks (SFN) [|. The phenomenon of Percola- 
tion has also been studied on such networks 0, 0] . In ad- 
dition disease spreading models like susceptible-infected- 
susceptible (SIS) 6] have been studied on networks for 
spreading diseases in the society or spreading of viruses 
in the Internet. 

Over the last few years it has been observed that the 
nodal degree distributions of many real- world networks, 
e.g., World Wide Web 7] and the Internet are char- 
acterized by power law tails: P{k) ^ k^^ (degree k of 
a node being the number of links attached t o it) . These 
networks are called 'scale-free networks' |^ llOl Hit ll^ 
due to the absence of a characteristic value for nodal de- 
grees. Theoretically a number of graphs are generated 
to model SFNs. One of them is by Barabasi and Albert 
(BA) which has the following ingredients, namely: (i) A 
network grows from an initial set of rrio nodes with links 
connecting all rrio nodes. At every time step a new node 
is introduced and is randomly linked to m(< nio) distinct 
previous nodes, (ii) Any of these m links introduced at 
time t connects a previous node i with an attachment 
probability tt^ (t) which is linearly proportional to the de- 
gree ki{t) of the j-th node at time t: Hi{t) oc ki{t). For 
BA model 7 = 30. 

Sandpile models are the prototype models of Self- 
organized Criticality (SOC) [ll El Hi ES E3 • In these 
models, long ranged correlations both in space as well 
as in time spontaneously emerge under a self-organizing 
dynamics, in absence of a fine tuning parameter. In its 
very general form, a sandpile model can be defined on 
an arbitrary connected graph, having a set of vertices 
connected by another set of edges. An integer height 
variable hi representing the number of grains in the sand 
column is associated with every vertex i of the graph. 
Starting from an arbitrary initial sand height distribu- 
tion the system is driven by adding unit grains of sand 
at the randomly selected vertices hi hi + 1. This 
sand column is said to be unstable when the height hi 
exceeds a pre- assigned threshold value he- An unstable 



sand column must topple and in a toppling it looses some 
grains which are distributed among the neighboring sites 
'l3j. This creates an avalanche of sand column topplings 
and the strength of such activity measures the size of the 
avalanche. There must be some 'sinks', i.e., a set of ver- 
tices through which grains flow out of the system so that 
in the steady state the balance of fluxes of inflow and 
outflow currents is maintained. 

In this paper, we have studied the deterministic sand- 
pile model on a scale-free network placed on an Euclidean 
substrate, namely a square lattice. The motivation of this 
study is to acquire support for the validity of our recent 
conjecture 0| that in a sandpile model the precise bal- 
ance at all lattice sites (except on the boundary) between 
the number of outflowing grains Hi which are distributed 
among the neighbouring sites in a toppling at the site i 
and the number of inflowing grains H'^ received by the 
site i when its all neighbouring sites topple once ensures 
that the sandpile model behaves like the BTW model 
|l3l | with a multiscaling avalanche size distribution. The 
absence of the site-to-site balance of Hi = H[ leads to 
the behaviour of Manna sandpile El ■ Below we define 
and study the sandpile model on the SFN where Hi is 
equal to the degree ki of the SFN and therefore is an ex- 
tremely fluctuating quantity. In spite of that the equality 
Hi — H[ is maintained by construction. We see below 
that the sandpile model on the optimized SFN indeed 
behaves like the BTW model. 

Recently, BTW sandpile model has been studied on a 
static model of SFN IQj . In contrast to the usual sand- 
pile models there are no specific sinks at fixed positions. 
Instead, during a toppling any grain can evaporate from 
the system from any arbitrary node with a small prob- 
ability /. The distribution of avalanche sizes (s) which 
do not dissipate (i.e., grains do not evaporate in these 
avalanches) is: 

Prob(s) ^ s""^ exp(— s/sc) (1) 

where the cut-off of the avalanche size Sc ~ !//• It is 
to be noted that the cut-off size does not depend on the 
network size N but only on the dissipation rate. 

To claim that a dynamical process active in a system is 
self-organized critical, it is important to ensure that both 
long ranged spatial and temporal correlations dynami- 
cally evolve in this system. For the ordinary BTW or 
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FIG. 1: Scaling of the probability distribution of avalanche 
sizes for the BTW model on the un-optimized SFN for differ- 
ent system sizes: L = 256, 512 and 1024. The data collapse 
gives the values of the scaling exponents D = 2 and /? = 3.2 
giving the avalanche size exponent '^un—opt — l.D. 

Manna sandpilcs on systems of spatial extension L this is 
verified in the following ways: (i) The avalanche size dis- 
tribution has a power law distribution Prob(s, L) ~ 
for some intermediate range and this range should in- 
crease with the system size as the cut-off of the avalanche 
size distribution increases as Sc{L) ~ L^. (ii) The aver- 
age size of the avalanches increases with the system size 
L, {s{L)) ~ L" and ly = 2 since the grains while exe- 
cuting a diffusive motion have to travel distances of the 
order of L to go out of the system through sinks situated 
on the boundary Non-zero values of v and D indicate 
that system has avalanches of all length scales and the 
process is indeed critical. 

If a sandpile is grown in a closed system (i.e., a system 
which has no sinks and grains do not evaporate from this 
system) the system eventually reaches a state when an 
"infinite avalanche" which continues for ever and never 
stops. Now, if a slow dissipation rate is introduced like 
every 1// topplings one grain is dissipated from any ar- 
bitrary site of the system, there is no infinite avalanche, 
the system indeed reaches a stationary state, but the 
avalanche sizes arc no more of all length scales. This 
is because the large avalanches loose their strengths by 
dissipation of grains. 

For a network or a graph in general, there is no con- 
cept of space, only the connections by links between the 
nodes. One still can define a distance between an arbi- 
trary pair of nodes on a network measured by the number 
of links on the shortest path connecting the two points. 
The largest of all possible shortest paths is called the 
diameter of the network. A small world network has 
the diameter varying logarithmically with the number 
of nodes: V{N) oc log A?". Since scale-free networks are 
small world networks, it is difficult to observe long ranged 
spatial correlations in sandpile model on SFNs. 




FIG. 2: A cost optimized Barabasi- Albert scale- free network 
on a square lattice of size L — 64. The SFN of A = 
nodes is generated by usual BA algorithm whose nodes are 
randomly assigned lattice sites. The cost function C (the total 
wiring length) is then minimized by a large number of trials 
as described in the text keeping the nodal degree distribution 
intact. Large degree nodes are visible. 



Here we study the sandpile model on a SFN con- 
structed on a square lattice of size L x L. We first con- 
struct a BA SFN oi N = L'^ nodes. The network starts 
growing with an initial set of TOq = {fn + 1) nodes. Each 
of these nodes is linked to all other m nodes forming a 
(m -|- 1)— clique. After that new nodes are added to the 
network one by one and each suc:li node is connected to 
m randomly selected distinct nodes of the already grown 
network with probability tt^ (t) . This process stops when 
the network size has grown to nodes. In our calcula- 
tion we use m = 2, therefore our network has nodes, 
2L^ — 3 links among the nodes and has many loops. The 
nodes of the network are then assigned randomly with 
uniform probability the sites of the square lattice. If two 
nodes are linked, the corresponding lattice sites are con- 
nected by straight lines. Thus we place the BA SFN on 
the square lattice. 

Clearly the degree distribution of such an Euclidean 
SFN is exactly the same as that of the BA SFN. To study 
the sandpile model we assume that each site (except for 
sites on the boundary) has a site dependent critical height 
/i? of stability which is equal to the degree ki of the node 
at that site. Therefore in a toppling, the sand height at 
site i is reduced to: hi ^ hi — hi and in a deterministic 
toppling dynamics like BTW model, all the fcj neighbours 
receive one grain each. The outlet of the system is at the 
boundary. Therefore every boundary site (except the cor- 
ner sites) has the threshold heights /i? = fc, -|- 1. This 
implies that in a toppling at the boundary site one grain 
goes out of the system and never comes back. Similarly 
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FIG. 3: Scaling of data for the probability distribution of 
avalanche sizes for the BTW model on optimized SFN and 
for different system sizes: L = 32, 64 and 128. The data 
collapse gives the values of the scaling exponents D = 2 and 
P — 2A giving the avalanche size exponent Topt ~ 1.2. 

at the corner sites the threshold heights are h'i — ki + 2. 
Such mechanism of outflow of grains through the bound- 
ary sites guarantees that the sandpile dynamics on the 
Euclidean SFN must reach a stationary state. 

In this Euclidean SFN any site is connected to any 
other site with equal probability and therefore the av- 
erage link length {£ij) is large and of the order of the 
system size L. In a toppling the grains therefore jump 
large distances on the average. We first study a determin- 
istic sandpile model on such a network. In this sandpile 
model a grain jumps a distance around L in a toppling. 
Consequently, the spatial extent of all avalanches, small 
or big, are around L. 

The avalanche size is measured in two ways: (i) the 
total number of topplings s and the number of distinct 
sites toppled a. A power law distribution of the avalanche 
sizes with a simple finite size scaling implies that the dis- 
tribution function Prob(s, L) obeys the following scaling 
form: 

Proh{s,L)^L-Pf{^), (2) 

where the scaling function f{x) ^ x^'^ in the limit of 
X giving r = r and D are the exponents 

of the avalanche size distribution. One immediate way 
to check validity of Eqn. (1) is to attempt a data col- 
lapse by plotting L'^Prob vs. s/L^ with trial values of 
the exponents. It is now well known in the literature 
that for the Manna stochastic sandpile model the distri- 
bution obeys ESS with TManna ~ 1-28 where as for the 
BTW sandpile the probability distribution Prob(s, L) of 
this measure has been found recently to obey a multi- 
scaling ansatz . In Eig. 1 we plot the scaling of 
the avalanche size distribution data for the BTW sand- 
pile on the un-optimized BA SEN on the square lattice 
for L = 256, 512 and 1024. The best collapse works 
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FIG. 4: Comparison of the BTW model on ordinary square 
lattice (solid line) and the BTW model on the optimized SFN 
on square lattice (dotted line): (a) the moment exponents 
(7(g) vs. q and (b) da{q)/dq vs. q. 

for D = 2 and /? = 3.2, giving Tun-opt = [i/D = 1.6. 
The stochastic Manna sandpile is also studied on the un- 
optimized BA SEN on the square lattice for system sizes 
upto L=1024 again. We estimated D = 2 and f3 = 3.0 
giving Tuuapt = 1-5. We believe that Tun-opt for both 
the BTW and Manna sandpiles on the un-optimized BA 
SFNs are indeed mean-field like and both the exponents 
should be actually 1.5. Similar slight deviation from 1.5 
was also observed in for the deterministic case. It is 
also observed that the area a has a similar distribution. 

Recently a cost optimized SEN on Euclidean space has 
been constructed [20j where the total sum of the link 
lengths is optimized keeping the nodal degree distribu- 
tion exactly same as that of the original SEN. For such 
a construction one defines a cost function C{N) as the 
total wiring length in terms of the symmetric adjacency 
matrix A of size N x N (which has elements = 1 if 
there is a link between the pair of nodes i and j and 
otherwise) and the distance £ij between nodes i and j 
as C{N) — Hiyja.ijlij. The optimization process is es- 
sentially a rewiring process maintaining the nodal degree 
distribution intact. It starts with a BA SFN constructed 
on a square lattice as mentioned above. A pair of distinct 
links of the SFN is chosen whose nodes are not linked oth- 
erwise. One end of each link is then opened and rewired 
suitably to another node of the quartet so that total sum 
of the rewired length is smaller. More precisely, the first 
node Til is randomly selected from the set of N nodes 
and the second node n2 is randomly selected from the ki 
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neighbours of ni. In the same way 7^3 ni 7^ n.2) is se- 
lected randomly from N nodes and n^l^ ni ^ ni) is cho- 
sen from fca neighbours of ria. Clearly this move conserves 
the link numbers as well as degree distribution. Rewiring 
is done following this decision: If both nina and 712/14 are 
not linked and also £12 + £34 is greater than ^13 -I- £24 we 
link n\n-i and n2n4. Another possibility is if nin4 and 
7^27^3 are not hnked but t\2 + ^34 is greater than £14 -I- ^23 
then we link 711/14 and 712/13. If both cases are possible 
we accept one of them with probability 1/2. If only one 
is satisfied we accept that. After rewiring we remove the 
links 711/12 and 712714. If none of the two is satisfied we go 
for a fresh trial. On repeated trials of these moves the 
cost function gradually decreases. Initially it decreases 
very fast but eventually the success rate becomes very 
slow. To monitor the optimization process we kept track 
of the average link length. Our best possible effort yields 
the average link length {lij) ~ 1.75 lattice constant. A 
picture of the optimized network is given in Fig. 2. In 
this best possible optimized network the link lenghs I 
have an exponential distribution as: Vit) ~ exp(— g£) 
with g Ri 1.16. Also the diameter of the network 'D{N) 
is measured and is observed to grow as where fi is 
estimated to be 0.40 ± 0.02. Therefore this network is 
scale-free but not a small-world network. 

The deterministic BTW sandpile model is then studied 
on such a network. The avalanche size distribution is 
calculated for three different system sizes L = 32, 64 and 
128. It was difficult to go beyond this size because of 
the large optimization times required. First we tried to 
make a scaling plot of the size distribution data. In Fig. 3 
we show this plot, which shows reasonably well collapse 
of the data in the intermediate range of the avalanche 
sizes. The corresponding /3 and D values fitted are 2.4 
and 2 respectively giving a possible value of Topt ~ 1.2. 
However for large avalanche sizes the collapse is much 
worse and the data for different system sizes separate 
out. This is a typical behaviour of the BTW like models 
which show strong presence of the multiscaling behaviour 

mill. 

The mutiscaling behaviour is studied in more detail by 
the evaluation of the various moments of the avalanche 
size probability distribution. The q-th moment of the 
distribution is defined as (s*) = / s'^Prob(s, L)cis. Incase 
the distribution Prob(s, L) obeys the finite size scaling 
behaviour for the whole range of avalanche sizes, it can 
be shown that (s"^) ~ L'^'-i^ where a{q) = D{q — T + 1) for 
q > r — 1 and cr{q) = for < q < r— 1. The q dependent 
exponent <T{q) is determined from the slope of the plot of 
log(s''(L)) with logL, which in our case are for L = 32, 64 
and 128. The interval between successive q values is 0.02 
and moments are calculated at 251 values of q between 
and 5. In Fig. 4(a) we show a plot of cr(g) vs. g on a linear 
scale. In Fig. 4(b) the derivative of a{q) is plotted with 
q. Had the Prob(s, L) followed a simple FSS behaviour 
the da{q)/dq in Fig. 4(b) would have saturated for large 



q values. In stead, the curve gradually increases with 
q, very similar to the multiscaling behaviour of BTW 
model. To compare we plot both a{q) and da{q) / dq of the 
ordinary BTW on square lattice studied for same system 
sizes with different line styles. We see that in both plots 
the behaviour is very similar and the difference between 
the two curves is very small, within 2-3 %. 

The stochastic Manna sandpile is also studied on the 
optimized SFNs for small system sizes L 32, 64 and 128. 
The scaling exponents are estimated as: D — 2.62 and 
/3 = 3.4 giving Topt ~ 1-3. This value of t is compared 
with the corresponding r « 1.28 value of the ordinary 
Manna sandpile. 

To summarize, we studied the BTW sandpile model 
on a Barabasi- Albert scale-free network oi N = nodes 
where the nodes are the sites of a square lattice of size L. 
The SFN is then optimized minimizing the total wiring 
length but keeping the degree distribution intact. On 
such an optimized SFN on the Euclidean space we ob- 
serve that the sandpile model has the same scaling be- 
haviour as the BTW model where as the deterministic 
sandpile on the un-optimized SFN has a mean-field like 
behaviour. 
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